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In the intermediate state of a thin type-I superconductor magnetic flux penetrates in a disordered 
set of highly branched and fingered macroscopic domains. To understand these shapes, we study in 
detail a recently proposed "current-loop" (CL) model [R.E. Goldstein, D.P. Jackson, A.T. Dorsey, 
Phys. Rev. Lett. 76, 3818 (1996)] that models the intermediate state as a collection of tense current 
ribbons flowing along the superconducting-normal interfaces and subject to the constraint of global 
flux conservation. The validity of this model is tested through a detailed reanalysis of Landau's 
original conformal mapping treatment of the laminar state, in which the superconductor-normal 
interfaces are flared within the slab, and of a closely-related straight-lamina model. A simplified 
dynamical model is described that elucidates the nature of possible shape instabilities of flux stripes 
and stripe arrays, and numerical studies of the highly nonlinear regime of those instabilities demon- 
strate patterns like those seen experimentally. Of particular interest is the buckling instability 
commonly seen in the intermediate state. The free-boundary approach further allows for a calcu- 
lation of the elastic properties of the laminar state, which closely resembles that of smectic liquid 
crystals. We suggest several new experiments to explore of flux domain shape instabilities, including 
an Eckhaus instability induced by changing the out-of-plane magnetic field, and an analog of the 
Helfrich-Hurault instability of smectics induced by an in-plane field. 



I. INTRODUCTION 

A longstanding problem in macroscopic superconduc- 
tivity is that of understanding the complex patterns of 
flux penetration observed in the intermediate state of a 
type-I superconductor. This state is observed when a 
thin superconducting slab is placed in a perpendicular 
magnetic field. Unlike type-II superconductors, where 
the field penetration is in the form of tubes each with a 
quantum of magnetic flux, type-I systems are observed to 
form intricately branched and fingered macroscopic flux 
domains Thus, instead of establishing a Meissner 

phase, in which the magnetic induction B = uniformly, 
the demagnetizing effects of the large aspect ratio force 
the sample to break up into regions, some of which are 
uniformly superconducting (with B = inside) and oth- 
ers that are normal (B ^ 0). Figure [l] shows a typical 
example of these patterns |3|. The superconducting re- 
gions appear black, having been decorated with a pow- 
der (niobium) that is itself superconducting at the sample 
temperature and thus migrates to the regions of low mag- 
netic field. Other imaging techniques include Hall probes 
H and magneto-optics 0] . All reveal similar structures. 

The sample in Fig. Wis at an applied magnetic field 
H a that is very close to the critical field H c at which 



the sample would be completely normal, so the minority 
phase is superconducting. Similar patterns are observed 
when H a / H c is very small, but now the minority phase 




FIG. 1. The intermediate state of a thin slab of indium, in 
which the superconducting regions (black) are decorated with 
niobium (black). The applied field H a is close to the critical 
field H c (h = H a /H c = 0.931). Adapted from Haenssler and 
Rindcrcr ||. 



1 



is normal; the sample consists of fingered and branched 
flux domains in a matrix of superconductor. These do- 
mains have a characteristic field-dependent finger width, 
and the branched domains have three-fold vertices. 

For many years it has been known || that more or- 
dered flux domain states may be observed in the pres- 
ence of a small in-plane component of the applied mag- 
netic field. It is also known that the domain morphology 
is not a thermodynamic state function; it depends on 
the path in field-temperature space through which the 
sample has been brought to a given point 0. Thus, for 
instance, cooling in zero field below the transition tem- 
perature and then applying the field tends to produce 
patterns in which normal domains are embedded in a 
matrix of superconductor, whereas when the same point 
in T — H space is reached by cooling below T c in a fixed 
field the normal domains connect to the sample edges jp . 
These observations suggest that the patterns are not true 
ground states of the system — the sample is not in true 
thermodynamic equilibrium. 

Despite the ubiquity of these patterns, there has until 
recently been no theoretical explanation for their form. 
The earliest attempt, prior to the detailed experiments 
described above, was by Landau R], who considered the 
laminar state: a periodic array of alternating supercon- 
ducting and normal domains (Fig. ||). Exploiting the 
translational invariance of the pattern in the direction 
parallel to the stripes, the cross-sectional shape of the 
domain walls and the associated bending of the mag- 
netic field lines become purely two-dimensional problems 
amenable to conformal mapping techniques. The pri- 
mary result of this calculation is a determination of the 
laminar state free energy as a function of lamina spacing, 
applied field, superconductor-normal surface energy, and 
slab thickness. 




FIG. 2. The laminar state in a thin slab. The applied field 
H a is normal to the slab. An ordered laminar structure is 
obtained with an additional in-plane component Hi. 



The free energy so obtained depends on two param- 
eters: the area fraction occupied by the normal state 
and the repeat distance of the structure. For thick slabs, 
one finds to a good approximation that the equilibrium 
area fraction is set by the reduced field h — H a /H c . 



Deep within the slab the domain walls are aligned with 
the field, but they flare along with the magnetic field 
lines near the surface. The equilibrium field-dependent 
width has been found to be in good agreement with sub- 
sequent experiments Stability calculations about 
this state are, however, precluded by the reliance on con- 
formal mapping techniques, and no systematic calcula- 
tions of this type have been performed. Subsequently, 
and in light of experiments showing nearly circular flux 
domains, there have been several calculations of the ener- 
gies of periodic arrays of simple geometric structures [|j . 
These studies generally find that the free energies of these 
periodic structures are very close to each other, so there 
is little thermodynamic driving force preferring one over 
the other. None of these theories offers an explanation 
for the appearance of branched flux domains. 

The dual features of (i) disordered domain shapes and 
(ii) path dependence to the patterns makes it clear that 
an understanding of the intermediate state requires a 
theoretical approach that (a) treats the superconductor- 
normal (SN) interfaces as a free-boundary problem, and 
(b) addresses the dynamics of those interfaces. That one 
should consider an interfacial representation at all, rather 
than a more microscopic approach based on equations of 
motion for the superconducting order parameter and the 
vector potential, is made clear by the strong separation 
of scales that exists between the domain size (typically 
fractions of a millimeter) and the width of the SN inter- 
face. For a strongly type-I superconductor the width is 
set by the coherence length £, which is on the order of 
0.2 jum. 

Several recent studies |9|,[l0| have emphasized strong 
connections between the motion of SN interfaces in the 
presence of a magnetic field and the dynamics of solid- 
liquid interfaces in the process of solidification. The key 
to this relationship is the fact that the magnetic field in 
the normal state obeys a diffusion equation JnJ, analo- 
gous to the diffusion of latent heat in the solidification 
problem. These suggest that a diffusional instability like 
the Mullins-Sekerka instability |^2| should occur in the 
process of flux invasion. Numerical studies of the time- 
dependent Ginzburg-Landau model confirmed these in- 
stabilities, which can lead to highly ramified domain pat- 
terns. 

It was also found by asymptotic methods |l3f| that in 
the limit of a sharp interface it is possible to integrate out 
the magnetic field contributions and arrive at a nonlocal 
free-boundary theory for the SN interface alone. The 
nonlocality is both temporal and spatial, the latter tak- 
ing the form of a Biot-Savart interaction between distant 
segments of the interface. The appearance of this type 
of coupling reminds us that the supercurrents that pro- 
vide the screening of the applied magnetic field flow along 
the SN interfaces. Many years ago Pearl H], and later 
Fetter and Hohenberg jlq] showed that the interactions 
between currents flowing in thin slabs are long-ranged. 
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Thus, vortices in a thin film interact with a potential 
whose long-range form is unscreened, V(r) *~ 1/r, while 
at short distances V(r) ~ ln(A/r), where A is a cutoff. 
Screening is unimportant at long distances because the 
interaction energy is dominated by the electromagnetic 
fields (in vacuo) above and below the slab. 

Based on all of these observations, we recently pro- 
posed [jl6| a "current-loop" (CL) model as an approxi- 
mate description of the intermediate state. The model 
describes the patterns as a collection of tense ribbons 
of current, interacting with the long-range Biot-Savart 
interaction of currents in free space, and subject to a 
constraint of global flux conservation. It is based on sev- 
eral simplifying assumptions. First, the domain walls 
are taken to be vertical, rather than flared as in Lan- 
dau's more precise calculation. Second, only the super- 
currents in those walls are considered — surface supercur- 
rents on the top and bottom of the sample are ignored. 
These approximations allow for an averaging process over 
the slab thickness that maps the simplified model onto 
one of self-interacting contours in the plane. Such free- 
boundary approaches are well-known for hydrodynamic 
problems such as Saffman- Taylor fingering llTj], as well 
as in systems governed by partial differential equations 
of the reaction-diffusion type. In the latter case, the dy- 
namics may be reduced to that of an interface when the 
scale of the patterns is large compared to the width of 
the front connecting two locally stable states. 

While several important phenomena are precluded 
from study in this model (such as domain fission), it has 
the virtue of providing a simple picture of the mechanism 
of shape instabilities in this system, and may form a use- 
ful starting point for a more precise treatment jlSfl . It 
also is strikingly similar to models for the energetics and 
dynamics of interfacial pattern formation in a variety of 
other systems that display "labyrinthine" patterns from 
a competition between interfacial tension and long-range 
electromagnetic interactions jl9) . These include fingering 
instabilities of magnetic fluids in Hele-Shaw flow po| P4| , 
Langmuir monolayers at the air-water interface |25| |27f| 
and thin magnetic films pq| . It is also equivalent to a 
reaction-diffusion system p9|-|32]| in which chemical fronts 
between two locally stable states may form space-filling 
disordered labyrinths similar to those seen in experiments 
on chemical pattern formation in a gel reactor [ ^3|j34| . 

In Section II we set the stage for a macroscopic 
model of the intermediate state by reviewing both the 
sharp-interface limit of the time-dependent Ginzburg- 
Landau model for nonequilibrium superconductivity and 
the derivation of long-range interactions between cur- 
rents in the slab geometry. The conformal mapping solu- 
tion to Landau's model of the intermediate state and the 
simpler straight-lamina model are derived in Section III, 
in which careful attention is paid to the consequences 
of those long-range interactions. The structure of the 
free energy as a function of stripe periodicity and width, 



surface tension and applied field in the two models is 
shown to be essentially equivalent, as are their predic- 
tions for the equilibrium stripe width as a function of 
field. The current-loop model is described and applied to 
the energetics of the laminar state in Section IV. Insta- 
bilities of single flux stripes and the elastic properties of 
the laminar state are found in Section V. The concluding 
Section VI summarizes the new experimental predictions 
that arise from the correspondence between the interme- 
diate state and others such as Langmuir monolayers and 
smectic liquid crystals. An appendix gathers together 
calculational details of the stability analyses and elastic- 
ity calculations. 

II. PROPERTIES OF THE MACROSCOPIC 
MODEL 

A. The sharp-interface limit 

Here we review the results of asymptotic methods ap- 
plied to the time-dependent Ginzburg-Landau (TDGL) 
model of superconductivity in the limit of sharp inter- 
faces. The microscopic parameters that enter the TDGL 
include the charge e* and mass m of a Cooper pair, a di- 
mensionless order parameter relaxation time 7, and the 
conductivity a of the normal phase. The coupled equa- 
tions of motion for the order parameter ip and the scalar 
and magnetic vector potentials, <fr and A, are 

+a-0 - b\i/j\ 2 ip, 
V x V x A = 4f (J„ + J s ) . (2.2) 

where J„ and J s are the normal and supercurrents, 

J n = o {-V</> - d t A) 

J. = r^(fW-W)-- M 2 A. (2.3) 

The bifurcation parameter a contains the important tem- 
perature dependence, a — ao (1 — T/T c ) and controls 

the correlation length £ = h/(2m\a\) 1 ^ 2 and penetra- 

1/2 

tion depth A = [m6/47re* 2 |a|] . Finally, the Ginzburg- 
Landau parameter k = A/£. 

When £ is small the superconductor-normal interface 
will be sharp, and a natural measure of smallness is pro- 
vided by the inverse distance from the critical point. Set- 
ting a = a/e 2 and rescaling all distances and fields with 
appropriate powers of e, one finds that both A and £ scale 
identically. The sharp-interface limit is then achieved by 
£, A — > with k fixed. It is thus possible in this limit to 
continue to distinguish between type-I (k < l/y/2) and 
type-II (k > l/\/2) systems. 
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Existing derivations of interface equations of motion 
in this context presume translational invariance in the 
direction along which the magnetic field is applied, and 
hence describe "bulk" superconductors JTjS. The results 
may be summarized as follows. Far away from the inter- 
face (the "outer" solution in the sense of matched asymp- 
totics) , the magnetic field in the normal region obeys the 
diffusion equation 



h t = DV 2 h , 



(2.4) 



with diffusion constant D = l/4na. The boundary con- 
dition on the magnitude hi of h on the interface C is 



h 



Jj2 [vSnK. 



(2.5) 



where <jsn is the interfacial tension (with dimensions of 
energy/area), K. is the interface curvature, T is a known 
kinetic coefficient, and v n = n • r t is the normal compo- 
nent of the interface velocity. 

The equation of motion is an expression for v n in terms 
of the function r(s, t), and it follows from a solution to 
the di ffusi on equation (|]J) given the boundary condi- 
tions (|2.5| ). As with all such problems, this involves a 
convolution over previous times and all space. But in 
the limit D — > oo, or vanishing normal state conductiv- 
ity, the temporal nonlocality disappears. This yields the 
equation of motion 



' r t {s,t) 



U 2 - IT 2 



8tt 

Hi 
8tt 2 



ds 



,R x t(s') 
~~R 2 



(2.6) 



with R = r(s) — r(s'). The dynamics ( |2.6| ) has the vari- 
ational form 



, 1 fflcff 

V9 Sr 



(2.7) 



where g is the interface metric and the effective interface 
Hamiltonian for a single domain is 

H 2 — H 2 
H cS [r] = c — -A + a SN L 

07T 



jds jds'i{s) ■ i(s')g(R) , (2. 



with G(R) = —(l/2n) InR is the Green's function of the 
Laplacian in two dimensions, and A and L the area and 
perimeter of the domain. As anticipated in the intro- 
duction, we see that the energy associated with the dis- 
tortion of the magnetic field lines threading a domain is 
represented by the self-induction of the boundary. The 
other contributions are simply the line energy <jsnL and 
an area term associated with the magnetic pressure and 
condensation energy. 



The limits of a sharp interface and zero normal state 
conductivity that produce a temporally local but spa- 
tially nonlocal contour dynamics have parallels in a sim- 
pler reaction-diffusion system of recent interest [^9] [32| . 
This is the FitzHugh-Nagumo model ]35| of the coupled 
dynamics of an activator u and inhibitor v, considered in 
two spatial dimensions, 



lit = DVu 
ev t 



F'{u)~ p(v-u) 



V 2 v 



(2.9) 



These partial differential equations are written in a 
rescaled form in which the activator diffusion constant 
D is normalized to that of the inhibitor, while e is a ratio 
of their characteristic times. The function F(u), whose 
derivative F'(u) appears above, is a double- well potential 
that describes the auto-catalytic behavior (and bistabil- 
ity) of the activator. When e ^ 0, the coupled dynamics 
is not a gradient flow in any standard form. 

It is clear that the two systems ( |2.2| ) and ( |2.9| ) share 
many features. In each case there is an order parameter 
field (ip or it) with an underlying bifurcation (the former 
continuous, the latter of first order), which is coupled to 
a diffusing field (unscreened or screened). In the case of 
the superconductor, the second field could be integrated 
out of the problem in exchange for an instantaneous non- 
local coupling of the field ip. An identical feature appears 
in the reaction-diffusion problem when the parameter e is 
small, for then the inhibitor relaxes on times short com- 
pared to that of the activator and will be slaved to u. 
Setting evt = 0, the field v may then be expressed as 
v = J Qu, where Q is the appropriate Green's function 
for the modified Helmholtz operator. This produces the 
nonlocal dynamics for u that is variational, u t = —8£/Su, 
with 

£[u} = Jd 2 x^D\Vu\ 2 + F{u)- l -p 2 

+ [ d 2 x [ dVu(x)g(|x - x'|)u(x') . (2.10) 



Now taking the limit D — > we may reduce (2.10) to 
a functional of the contours bounding regions in which 
the activator takes on values corresponding to different 
minima of F . That new functional (of the position vec- 
tors Ti(s) for each domain boundary, parameterized by 
arclength s) is @Q 



A£[{rJ]=7E^ + Ai71 E^ 



-y 

2 ^ 

i,3 



<j>ds j>ds' 



(2.11) 



with 7 the line tension, and where Li and Ai the perime- 
ter and enclosed area of each domain. The nonlocal term 
is again in the form of the self- and mutual-induction of 
(here fictitious) current loops encircling each of the do- 
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B. Long-range forces 



III. THE LAMINAR STATE 



As further motivation for the CL model, we recall 
Pearl's derivation Q of the long-range (unscreened) po- 
tential between point vortices in a thin film. The starting 
point is the solution for the supercurrent surrounding a 
single vortex located at the origin. In the slab geometry, 
the order parameter has the simple form 



V>(r) = Miry" (0<z<d) 



(2.12) 



where ipo is the far-field value of ip, which vanishes for 
z < and z > d. The function / (r) describes the vortex 
core structure, and we take / = 1 on the (large) scales 
of interest. The second of the Ginzburg-Landau pair of 
equations (|2.2|) is then 



V x V x A = 4ttJ s = -^r 



^0 

2nr 



8- A 



(2.13) 



If the film is sufficiently thin, then we may average over 
the film thickness d to obtain 



V x V x A = 



A 2 



•00 ( 
2itr 



5(z) 



(2.14) 



In this form, it is readily apparent from the appearance 
of the term — (dA/ X 2 )5(z) on the right-hand-side that 
the screening is confined solely to the slab. 

The solution to this can be obtained through the use 
of Fourier-Hankel transforms and yields the supercurrent 



A. Landau's free-boundary solution 

In this first part of our discussion, we calculate the 
shape and optimal spacing for an assumed laminar ge- 
ometry for the intermediate state. While our work re- 
produces Landau's original calculation , our method is 
quite different, and we can provide an explicit expression 
for the function f(h) (see below). We are including this 
material in our discussion since (1) Landau's derivation 
is, to our taste, a bit obscure, so we hope that the present 
derivation will clarify the techniques and inspire further 
work on extending the calculations to more complicated 
geometries; and (2) the "exact" result derived here, for 
flared normal lamina, can be compared against the re- 
sults of the current loop model discussed below, allowing 
us to "calibrate" the current loop model. 

We consider the simplest semi-infinite geometry here, 
which we have illustrated in Fig. |[ The surface of 
the material is in the x — y plane, occupying the region 
z < 0. The field is applied perpendicular to the sam- 
ple, so that H a = H a z. We will assume that there is an 
array of normal-superconducting laminae periodic in the 
x-direction, the laminae being straight in the ^-direction. 
The normal and superconducting laminae have asymp- 
totic widths a n and a s , respectively, so that the period 
of the entire structure is a = a s + a n . 



8ttA 2 



H 1 {r/A)-Y 1 {r/A)- 



(2.15) 



where Hi and Yj are Hankel functions, and the effective 
penetration depth for the thin film is 

A™. (2.16) 

The interaction potential of a second point vortex placed 
a distance r from the first is obtained by multiplying 
( |2. 15 ) by 4> and integrating, with the result 



V(r) 



8ttA 



[Ho (r/A) - Y (r/A)] 



(2.17) 



At long distances r/A ^ 1, this is an unscreened poten- 
tial, 



V{r) 



(2.18) 



47r 2 r ' 

whereas as r — *■ the familiar logarithmic interaction 
between vortices appears, 



V(r) 



Z^a^Ua 



e c r 



(2.19) 



with C being Euler's constant. We anticipate therefore 
that the interactions between Meissner currents in a thin 
type-I slab should have a similar long-range character, 
and now proceed to a detailed calculation. 




FIG. 3. Laminae in the intermediate state. Superconduct- 
ing regions are shaded. Reference points P n , Q n , P' n , etc. are 
discussed in text. 

We start by noting that the typical lamina spacing 
is large compared to the superconducting penetration 
depth and the coherence length, so that the interfaces 
separating the normal and superconducting phase are 
sharp. Then, following Landau, we can work with the 
macroscopic Maxwell equations. In the normal regions 
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we have VxB = 0,V-B = 0, which can be solved either 
by introducing the scalar potential 4> through B = — V0, 
so that V 2 = 0, or the vector potential A through 
B = V x A, which in the gauge V • A = satisfies 
V 2 A = 0. We thus need to solve Laplace's equation in 
the x — z plane. This immediately suggests the use of 
complex variable methods; indeed, for the laminar struc- 
ture the only nonzero component of the vector potential 
is A y , so we may introduce the complex potential 



w 



• + iif) = <f> + lAy 



so that the complex magnetic field is 

B = B, r — iB x = — — , C = x + iz. 



dw 



(3.1) 



(3.2) 



The boundary conditions are the following. First, as z — > 
oo the field becomes the uniform applied field H a z, so 
that 



(3.3) 



Since the magnetic field vanishes in the superconduct- 
ing regions, continuity of its normal component implies 
that the field is purely tangential both along the normal- 
superconducting interface P n Q n an d along the segment 
which we refer to as the "nose" (OP n ) (Fig. ||). In the 
first case, the assumption of local thermodynamic equi- 
librium at this interface implies that the magnitude H n of 
the tangential component is a constant (as yet unknown) , 
so that the vector potential A y is also constant along any 
one interface. Along the nose the magnetic field is par- 
allel to the surface (H = H x x) but with a magnitude 
that is no longer constant. The field H n along with the 
periods a n and a s are determined a posteriori by min- 
imizing the energy. However, these parameters are not 
completely independent, since by flux conservation we 
must have H a a — H n a n . Therefore for fixed external 
field H a the energy will be determined by a and a n . 



1. Exact determination of the lamina shape 

The position of the interface separating the normal 
and superconducting regions is not known a priori, and 
must be discovered in the process of solving the problem. 
Although this sounds like a formidable task, it is made 
easier by recognizing that our magnetostatics problem is 
formally equivalent to the flow of an ideal incompressible 
fluid around an array of plates, the plates being the noses 
of the laminae. The field lines would be the streamlines 
of the fluid; the normal-superconducting interfaces cor- 
respond to fluid streamlines which have separated from 
the flow behind the plate ( "free-streamlines" ) . This cor- 
respondence is outlined in more detail in Table |. The 
shape of the free-streamlines can be determined by us- 
ing the hodograph method pQiWlj . The idea is that while 



the field lines in the £-plane may be complicated, the 
representation of the field configuration in the w and H 
planes is quite simple; a suitable conformal transforma- 
tion which maps the w plane onto the H plane would 
then lead to a relation between w and dw/dC,, the solu- 
tion of which will determine the shape of the interface. 
Consider first the w-plane. Since the magnetic field is 
tangent to both the nose segment and to the normal- 
superconducting interface, the magnitude of the vector 
potential is constant on these segments, as well as on the 
centerline shown dashed in Fig. |^. Far from the sample 
(at a fixed z — > oo), we know that A y ~ H a x; therefore, 
for laminae separated by a distance a, the vector poten- 
tial on the interface of the n — th lamina is A y — H a na. 
The potential behaves as —zH a as z — ► oo, hence the 
lamina correspond to the positive half of the <f> axis. The 
field configuration in the w-plane is shown in Fig. |](a). 
Next we consider the -ff-plane; this is simpler if we intro- 
duce the normalized magnetic field 



n 



H 



1 dw 
H^'dC' 



(3.4) 



3aH, 



ait 



2aH, 



en* 



a) 



aH 



at* 



W 



Pi 



a 
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b) 



Q' 



Rea 



FIG. 4. The w- and cr-planes. 

On the interface we have |t;| = 1, and we define 8 to be 
the tangent angle (as in Fig. |J) as rj — —e~ ie and only 
the phase 6 changes as we traverse the interface from the 
point P n to Q n ; this produces a semi-circle in the ?7-plane. 
On the nose r] x = 0, so that the segment P' n P n maps 
onto a horizontal straight line in the r\ plane. Therefore, 
the region exterior to the superconducting laminae maps 
onto the interior of a semi-circle, as shown in Fig. 0(a) 



G 



(the hodograph for the field). 

We now need to map the w-plane onto the ?y-plane. 
This is most easily carried out by the following sequence 
of transformations. First, we take care of the periodicity 
in the w-plane by mapping it onto the cr-plane with the 
transformation 



_ 2ttw / H a a 

a = e 1 



1. 



(3.5) 



so that the cr-plane has a single cut along the positive 
real axis, as shown in Fig. [|(b). Next, map the cr-plane 
onto the lower half of the A-plane using 



A = 



1/2 



(3.6) 



as shown in Fig. ||(b). Here k is a constant to be deter- 
mined from the boundary conditions. The next step is to 
map the 77-plane onto the A-plane using the Joukowsky 
transformation 



A 



2 \ v 



(3.7) 



To determine k, we notice that as z - > oo, w — > — oo, so 
that a — — 1 and r\ = H a /H n . Introducing the notation 
h a = H a /H n , we then have 



2^2 



k 2 = 



(3.8) 



a) 



p> — ** Rer| 



b) 




FIG. 5. The r/- and A-planes. 



Combining Eqs. ( 3.5)-( 3.8), we have for the complex 
potential 



H a a 
2vr 



In 



(v 2 + h 2 a )( v 2 + K 2 ) 



(I + 77 2 ) 2 



(3.9) 



To determine the shape of the lamina, we use Eq. (3.4) 
to write 

d£ <i£ dw 1 1 dw(rf) 

dr\ dw dr\ H n rj dr\ 

a (1- h 2 a ) 2 1 

7T h n 



If 



(v 2 + ha 2 )(ri 2 + hl)(^ + l) 



(3.10) 



The integral can be performed by using a partial fraction 
expansion. The constant of integration is determined by 
recognizing that when rj — 0, z — 0. We finally have 
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This equation implicitly determines the magnetic field rj 
as a function of position £. The particular parameteriza- 
tion is different from Landau's |t]|38| , but can be shown to 
be equivalent. This result also appears to coincide with 
result of Fortini and Paumier (3^] , although the method 
of derivation is entirely different. The equivalent fluid 
problem would be the separated flow past a plate placed 
symmetrically in a channel; the solution to this problem 
is given in R ef. | f37| , p. 39, Eq. (25a), which is the same 
as Eq. (3.11) above. 

To calculate the lamina shape, recall that on the 



normal-superconducting interface rj = 



-i0 



with 9 



the tangent angle on the interface. Substituting into 
Eq. (3.11), and separating real and imaginary parts, af- 
ter a great deal of algebra we obtain the following two 
parametric equations for the position of the interface: 
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(3.12) 



Using these parametric equations we can calculate the 
lamina shapes for different applied fields h a \ some rep- 
resentative results are given in Fig. |. The width 2b of 
the nose can now be determined from ( |3.12 ) by setting 
9 = 0, with the result 



6=|(l-M-^(l-^)sin- 1 



2hg ' 



(3.13) 
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FIG. 6. Lamina shapes calculated from Eq. (3.12) for sev- 
eral values of the reduced field h a . 



2. Energy of the laminar structure 

The total free energy of the laminar state has con- 
tributions from the condensation energy, the magnetic 
field, and the SN interfaces. The condensation energy 
is just —H 2 /8ir times the volume occupied by the su- 
perconducting phase. If the normal-superconducting 
interfaces did not bend at all, this energy would be 
— (H 2 /8ir)(Na s )L y d, with N the total number of lam- 
inae, L y the length of the sample in the y-direction, and 
d the thickness of the sample [note that the sample has 
a total area A — (Na)L y ]. To this we must add the con- 
densation energy lost due to the thinning of the super- 
conducting regions near the sample surfaces. The con- 
densation energy per unit area is then 
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dz, (3.14) 



where the factor of 4 in front of the integral accounts for 
the 4 corners of the lamina. The integral may be written 
as an integral over the tangent angle 9: 



x(z) 



dz 
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a 2 h a (\ - hi) /■^ tan- 1 (pcos^)(l - cos 2 9) 



2ir 2 



cos 9(1+ p 2 cos 2 9) 



where the second line was obtained by substituting from 
( 3.12 ), and with p = 2/i a /(l— h 2 ). After a partial fraction 
ex pansio n and several integrations by parts the integral 
in ( 3.15| ) is found to be 



^ ^ln(p + q) - ^ In q 



(3.16) 



with q = (1 +P 2 ) 1 / 2 . We then have for the condensation 
energy 
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with 
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(3.18) 



Next, we need to find the magnetic field energy, ob- 
tained by integrating B 2 /8tt over all space, both inside 
and outside the sample. For our periodic structure the 
field energy per unit area becomes 
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(3.19) 



where the factor of 2 accounts for the top and bottom 
surfaces of the sample, and the integral is taken over 
the area C of one cell, as shown in Fig. [?]. We can write 
B 2 = {SI A y ) 2 in this two dimensional geometry, and then 



use the fact that V 2 A y = to write Eq. ( |3.19 ) as a line 
integral around the boundary dC of the unit cell: 

with A y (s) the vector potential on the boundary, and 
B s (s) the component of the magnetic field which is tan- 
gent to the boundary. 




FIG. 7. Contour for energy calculation. 

The advantage of this representation of the field en- 
ergy is that the vector potential is constant on the bound- 
aries, and it can therefore be brought outside the integral. 
Let's consider the various contributions to the integral. 
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On the normal-superconducting interface (the segments 
Q'P' and PQ) the integral J B s ds vanishes, since the 
field points in the same direction on the left half of the 
superconducting lamina as on the right half. On the mid- 
line between the laminae labeled n — 1 and n, the vector 
potential is H a a(n — l)/2 while on the next midline up it 
is H a (n+ l)/2; the integral J B s ds is equal in magnitude 
but opposite in sign for these two segments, since the in- 
tegration paths are in opposite directions. Adding these 
two contributions, and using B s = —dcf>/ds we have 



^ = ^ f B s ds 
A 4n J de 

= tk[<j>(a/2,-d/2) 



4>{a/2,L x /2)] 



(3.21) 



where L z is some large distance away from the top surface 
of the sample. The entire calculation of the field energy 
then reduces to finding the asymptotic behavior of the 
scalar potential along one of the streamlines (the midline, 
in this case) . By examining the behavior of Eq. (|3.9|) for 



the complex potential w, and Eq. (3.11) for the position 
£, as rj — > 1 and r\ — > h a , we find the asymptotic behavior 
<t>{a/2, L z /2) ~ -H a L z /2 + </>+, where 
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+2h a \nAh a . 



(3.22) 
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Substituting into ( 3.21 ), we have 



E m _ H 2 



H n H a H , . 

__ rf+ _ a/mag (/ la ) 



(3.23) 



(3.24) 



with 



/mag(^a) 



K 

2n 
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(3.25) 



The first term in Eq. (3.24) is the energy of the external 
field in the absence of the sample, which is of no inter- 
est and will be dropped from now on. The second term 
is the bulk magnetic field energy for a uniformly magne- 
tized sample pOfl , and the third term is the energy arising 



from demagnetizing fields (due to the partitioning of the 
sample into domains). 

Finally, we need to calculate the surface energy due 
to the normal-superconducting interfaces. If <jsn is the 
surface tension for the normal-superconducting interface, 
then the energy for a single interface is asNdL y [ pl| . 
Since there are two interfaces per lamina, and N lam- 
ina in the sample, the total energy due to the inter- 
faces is 2<jsN<lLyN = 2asNA{d/ a). We can introduce 
a length A, which is essentially the width of the inter- 
faces, through <7sn = (Hl/8ir)A; then the energy per 
unit area due to the interfaces is 



£int = Hi 2Ad 
A 8i o ' 



(3.26) 



Ad ding to gethe r all of the contributions to the energy, 
Eqs. (3.17), ( |3.24 ), and (3.26), and using the flux conser- 

H a a, we find 
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(3.27) 



with h = H a / H c . The energy must be minimized with re- 
spect to both a and a n . This results in very cumbersome 
expressions. Instead, we will minimize the first term in 
brackets with respect to a„, which yields a n — ha, so that 
H n = H c (and h a = h). This is the result used by Lan- 
dau, which is reasonably accurate as long as the surface 
and demagnetizing energies are small. If we substitute 
this back into the energy, we obtain 
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(3.29) 



This function is plotted in Fig. |8J. Its asymptotic be- 
havior as h — > is of interest in comparison with other 
approaches discussed below, and has the form 



h 2 , /0.56 

v ln hr 



(3.30) 



Finally, the equilibrium laminar period is obtained 
simply by minimizing with respect to a, yielding 
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(3.31) 
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FIG. 8. The functions fh(h) in Landau's model and f s (h) 
in the straight-lamina approximation. 

This implies that the spacing diverges for small h as 



ln(0.56//i) 



1/2 



(3.32) 



B. Energy of a straight-lamina model 

It is useful to compare the results from the Landau 
model to an alternative straight-lamina model, which is 
illustrated in Fig. ^|. In this model the field is still tangent 
to the normal-superconducting interface, but the magni- 
tude of the field is not constant along the interface. The 
magnetic field and the complex potential can be obtained 
using standard conformal mapping methods; this prob- 
lem is equivalent that of an ideal fluid flowing in a channel 
with an abrupt step pqj. The solution is 
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Using this solution it is possible to calculate the total 
energy, as in the Landau model. Making the simplifying 
assumption that H n = H c , we obtain 
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FIG. 9. Geometry of the straight-laminae model. 



with 



2tt 



fs(h) = 7- (1 + hf ln(l + h) - (1 - hf ln(l - h) 



-2h\n4h 



(3.36) 



This function is plotted in Fig. g for comparison with 
Landau's result (3.2S). At small h it behaves as 



h\ /0.68 



(3.37) 



which is very close to the Landau result (3.30). We see 



that the total energy for this model has the same qualita- 
tive dependence upon the lamina spacing a as the Landau 
model, although the Landau model has a lower energy for 
any reduced field h. 



C. Laminae shapes in finite-thickness plates 

In the analyses above we have assumed that the slab 
of superconductor is sufficiently thick that the shapes of 
the lamina walls can be computed as for a semi-infinite 
slab. When the thickness d becomes small enough the 
wall shapes will change. From the asymptotic behavior 



of the semi-infinite solutions (3.12), we deduce that the 
thickness approaches its asymptotic values for z — > — oo 



^(1 - h a ) - x{z) ~ exp(TTz/ah a ) . (3.38) 

The decay length ah a /n should then determine when 
finite-slab thickness effects become important. The 



asymptotic results (3.3C) and (3.37) show that the prod- 
uct ah vanishes very slowly (logarithmically) as h — > 0, 
so that while such finite-thickness effects become impor- 
tant in that limit, practically the relevant fields are very 
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small. From the asymptotics, the crossover field h x 
ah a /it = d is on the order of 



for 



h x ~ 0.56 exp(— A/nd) 



(3.39) 



For h < h x the slab thickness has no significant effect 
on the domain wall shapes. In the fluid dynamical anal- 
ogy, the finite-thickness slab calculation is equivalent to 
so-called Riabouchinsky flows around two plates, the de- 
tails of which the interested reader will find in standard 



references |37 



IV. THE CURRENT LOOP MODEL 

The analysis of the laminar state in the previous sec- 
tion has shown that in accounting for the flaring of the 
normal domains Landau's free-boundary approach yields 
a lower energy structure than a model with straight walls. 
But the analytical and numerical differences between the 
two approaches are relatively minor. In both models the 
supercurrents flow along the SN interfaces and on the top 
and bottom surfaces of the sample fll8| , [42"| ] . Just as the 
magnetic field in a solenoid is more nearly uniform when 
it is a tall thin cylinder than when it is short and wide, so 
too do we expect that the contributions from circulating 
currents along the SN interfaces will dominate when the 
flux domains are narrow and tall, at low applied fields 
h. This suggests that the basic physics of the laminar 
state can be understood by considering those circulating 
currents alone. We develop this current-loop model in 
the present section and show that it rather accurately 
reproduces the results of Landau's model. This calibra- 
tion is an important test of an approach that can easily 
be generalized to SN interfaces of arbitrary shape. 



A. Energetics of the current loop model 

For the purposes of this model, the intermediate state 
is described in macroscopic terms. As above, it has thick- 
ness d, total area A, and volume V = Ad, but now the 
SN interfaces encircling each of the normal regions be- 
long to a set {T>i} each with area A, and perimeter Lj. 
The two phases occupy volumes V s and V n — d^ Ji A il 
with V s + V n = V . Parameterizing each boundary by 
arclength, the position vectors of the interfaces are rj(s). 
As in the straight-lamina model, we assume that is in- 
dependent of z, neglecting the flaring of the domain walls 
near the film surfaces that is seen in Fig. ^. 

The total energy of the system is as before a sum of 
the condensation energy E c , the interfacial energy i?i n t, 
and the magnetic field energy E m , 



E[{r A] =E C + E i 



hit 



(4.1) 




FIG. 10. A current loop. 



And as before, 
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E c = V^Pn 
87T 



•Eint = 
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(4.2) 



where p n — A n /A is the area fraction for the normal 
phase, and where E c is measured with respect to the 
purely superconducting state. Global flux conservation 
relates the field in the normal regions to the applied 
field, H n = H a / p n , and by the tangential continuity of 
H across a S-N interface the field in the superconduct- 
ing region is H s = H„. The superconducting regions 
are perfectly diamagnetic; they each have a magnetiza- 
tion M = — (H n /ATr)e z , related in the usual way to the 
(Meissner) currents that flow along the SN (and top and 
bottom sample) boundaries. 

We compute the magnetic field energy as a sum of two 
contributions, the first of which is that of the domain 
magnetizations M in the presence of the external field 
(— J d 3 rH a • M). The second contribution is the self- 
and mutual-induction of the current ribbons. Expressing 
these in terms of the macroscopic quantities p n , etc., and 
the current-ribbon positions we have 



E n 



= y HaHn h 

47T V 



~\ M2 Y. J dz J dz 'j ds j ds ' 



Rii 



(4.3) 



where M — —H n /Air. Here, the vectors t; = tj(s) are 
unit tangents to the current ribbons and label the direc- 
tion of the current flow. By the usual screening processes 
in superconductors, the direction of the flow is so as to 
cancel the applied field in the superconducting regions 
and augment it in the normal regions (see Fig. |To| ) . The 
scalar product of the tangent vectors is however invari- 
ant under the overall reversal of the parameterizations 
(s — > — s). The current-current interaction is Coulombic, 
with Rij = {[r, t {s) - rj(s')} 2 + (z - z') 2 } 1 / 2 . While the 
z and z' integrals are readily performed (see below), the 
more elementary form (4.3) serves to remind us that the 
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current-current interactions are in their free-space form. 
By performing the z and z' integrals, the magnetic field 
energy becomes 



BL = V- 



.H„H„ 



-in 



(1 - Pn) 



where now the elementary free boundaries are contours 
in the plane, interacting with the potential 



$(R/d) = — [ dz [ dz' 
2« Jo Jo 



R 2 + {z- z'Y 



-1/2 



= smh" 1 {d/R) + R/d - V 1 + (R/d) 2 , (4.5) 

where R = |R| with R(s, s') = r(s) — r(s') the in-plane 
vector between points labeled by s and s' . As discussed 
elsewhere |23 this potential is Coulombic for R ^> d, 
$ d/(2R), but for R < d, $ w ln^e^d/fl), with the 
film thickness d acting as a cutoff. Note t he int eresting 
parallel with Pearl's interaction potential ( [2.17j ) among 
vortices in thin films. The fact that this interaction po- 
tential is in some cases identical and in others nearly 
identical to that found in the free-boundary approach to 
a number of other systems offers an explanation for the 
similarity in their behavior. Table II summarizes these 
analogies between the different systems. 
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FIG. 11. The laminar state as a collection of current loops. 
Arrows indicate directions of the current, tangential to the 
SN interfaces. 



B. Current-loop description of the laminar state 

Now we calculate the energy of the laminar state in the 
CL model. As before, the periodicity length in the plane 
is a, and the width of the normal lamina is a n ] we will 
assume that the pattern is periodic in the ^-direction. 
The nonlocal magnetic contribution is the only difficult 
one, and it proves more convenient to return to the orig- 
inal self-induction form of the tangent vector coupling, 
before averaging over the slab thickness. It is then easy 



to introduce a Fourier representation, which for a uni- 
form laminar structure yields the magnetic field energy 
E m per unit area 

a 
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Several straightforward integrations reduce this to 



(4.6) 
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aM 2 ^ [1 — cos {2irna n / a)] 
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x f^d +e -2,dn/a_ i y (4 7) 

Note that the last term contains all of the d-dependcncc 
of the sum. The leading contribution in the limit of large 
slab thickness is a "bulk" contribution expressible simply 
in terms of the stripe dimensions 



71=1 



[1 - cos (2ima n /a)) 2 a n (a-a„) , . 
Z2 = 77 12 • ( 4 - 8 ) 



It follows that the form of the free energy is exactly like 
that in Eqs. ( 3.27 ) and ( 3.2S ), but with a new function 
/, 



1 sin (mrh) 



1 



-2nnd/ a 



(4.9) 



As alluded to in the discussion of the laminar shapes 
in Landau's calculation, finite-thickness effects show up 
when d is comparable to a. 
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FIG. 12. Comparison of current loop model and Landau 
model for laminar state. Equilibrium laminar spacing a{h) 
for Landau model (solid) and current-loop model (dashed). 
Inset shows the functions /l (solid) and Jcl (dashed). 
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There are several noteworthy features of the function 
fch(h), particularly in the limit d/a ^> 1 which we con- 
sider here. First, as shown in Fig. [l^, it is rather close to 
the Landau function, and hence its implications for the 
equilibrium laminar thickness are in reasonable accord 
with experiment. Second, analytically it has the same 
structure as /l as h — ► 0, 



2tt 



In 



0.71 



(4.10) 



Third, it has an exact symmetry under the transforma- 
tion h — > 1 — h, a reflection in part of the straight SN 
interfaces presumed in the model. This symmetry is ab- 
sent in Landau's calculation and in the straight-lamina 
model by virtue of the currents on the slab faces. Fi- 



nally, the form of the magnetic field energy in Eq. (4.7) 
is identical to the field energy of a stripe array in a thin 
ferromagnetic film iflSlJi" 



C. A dynamical model 

The energetic competition between surface tension and 
self-induction present in the current-loop model appears 
in a number of other contexts (Table II), where it has 
been shown to produce also a rich dynamical behavior. 
While the precise connection between the Young-Laplace 
and Biot-Savart forces and the interface dynamics de- 
pends on the physical setting (e.g. Hele-Shaw flow with 
Darcy's Law, surfactant monolayers at the air-water in- 
terface with coupling to the fluid phase, reaction-diffusion 
systems), the phenomenon of repeated branching insta- 
bilities producing disordered lamellar structures is ubiq- 
uitous. This suggests that much can be learned by study- 
ing the very simplest dynamical law for interface motion, 
the local dissipation model p3[ , in which a local viscous 
drag acting at the interface balances the local pressure 
difference, the latter computed as a functional derivative 
as in Eq. (Oh. 



As a first step toward a full study of the many-interface 
current-loop model of the intermediate state, we study 
here the simplest mean-field description of a single cur- 
rent loop. That loop is assigned to a cell (analogous to a 
Wigner-Seitz cell) of area ^4 C eii from which we compute 
the area fraction p n — A n / A cc \\. In this approximate de- 
scription, the self-induction of the loop is retained in its 
full form, but the mutual induction term in the energy 
associated with all other loops only contributes a bulk en- 
ergy term like that seen in the laminar calculation (4.8). 
Moreover, the amplitude of the circulating currents is 
taken to be set by H c rather than the actual local field. 
This is equivalent to assuming that the actual area frac- 
tion departs little from its equilibrium value. The system 
is then characterized by a single geometrical parameter 
p = 2R cqu i v /d describing the aspect ratio, where i? oq uiv 



is the radius of the circle whose area is that of the initial 
condition, and a single energetic parameter, the reduced 
magnetic field h. All other parameters simply rescale 
time. 

Within this model, the normal component of the inter- 
face velocity is found by functional differentiation: 



n • r t (s) 
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n - A/C(s) 



-^^d S 'Rxt(,s')*(i?/d)| , (4.11) 

with JC(s) the curvature, and = $'(£) = 1 - (1 + 

£- 2 )!/ 2 j s the generalization of the Coulombic form of the 
Biot-Savart force to finite- thickness slabs (^3|. Finally, 
the pressure term is 



n = h 2 / P l - i 



(4.12) 



The kinetic coefficient rj may be estimated fl3|| from re- 
sults on the bulk properties of strongly type-I systems, 



8tt 8k B T c e ' 



(4.13) 



where again T c is the critical temperature and £o is t ne 
bare correlation length |4f| . 

A contour dynamics such as ( [1.11 ) is readily general- 
ized to account for surface tension anisotropy, a material 
feature that has long been suggested to play a role in the 
morphology of the intermediate state patterns Q, as it 
does in problems such as dendritic growth [^6|. When 
the anisotropy is g-fold, the parameter A has the form 



A = A [1 + e cos (q6)\ 



(4.14) 



Typical experiments show a(j = 4orq = 6 anisotropy 
p7[ . Our intuition suggests that the variation of o~sn 
through A will bias instabilities toward g-fold symmetry, 
and lead to preferred orientations of flux stripes produced 
from those instabilities. 



D. Instabilities; numerical studies 

Two regular geometries of flux domains have histori- 
cally been of interest: circles and stripes. In the next 
section we consider in detail the stability of stripes and 
stripe arrays (the laminar state) ; here we focus on finger- 
ing and branching instabilities of circular domains. Since 
linear stability analyses for circular interfaces have been 
presented in detail elsewhere in the context of closely re- 
lated models ||-|2|g||Q, we will not repeat them 
here in detail. Two important qualitative results from 
those studies are that for a given size (domain radius and 
slab thickness) (i) there exists a critical applied field be- 
low which the circle is stable and above which azimuthal 
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modes become active, and (ii) instabilities of increasing 
mode number occur with ever larger applied field. These 
properties may be illustrated through numerical studies 
of the contour dynamics, which allow us to see the highly 
nonlinear regime far beyond the instabilities. 




FIG. 13. Results from numerical studies of the current-loop 
model. Relaxation of a circular initial condition to a stable 
circular final state of different radius. Dashed line shows the 
area of the unit cell. Parameters are p — 5.0 and h — 0.1. 



An efficient numerical method for studying this non- 
local interface dynamics has been described in consid- 
erable detail elsewhere [ p7|]3(i|l . It uses pseudo-spectral 
techniques to solve for the time evolution of the local 
tangent angle 9{s), from which the (x(s),y(s)) coordi- 
nates of the interface are computed by basic differential 
geometry. For the purposes of verifying the analytical 
stability results as well as investigating such phenomena 
as mode competition, the initial contour is given a cur- 
vature K. perturbed from that of a circle, 

K(a) = — + [a n cos (na) + b n sin (na)] , (4.15) 

where Rq is the unperturbed radius and a — s/Rq. 

There are three basic phenomena that may be illus- 
trated with the contour dynamics. The first, shown in 
Fig. is the relaxation of a weakly perturbed circular 
domain whose initial area fraction is not the equilibrium 
value. This stable relaxation to a circle can occur if the 
applied field h is below the instability value at the as- 
pect ratio of interest. The asymptotic area fraction at 
long-times is p n ~ h, apart from a small correction due 
to surface tension. The approach of p n to this limiting 
value is shown in Fig. [14 
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FIG. 14. Time evolution of the normal area fracti on p „ 
for a single domain governed by the contour dynamics ( [f .ll|) . 
Curves (a), (b), and (c) correspond to Figs. [Dj, and |lq. 
Dashed lines indicate the relation p„ — h determined by the 
bulk energetic contributions alone. 

A second phenomenon occurs at higher h, and is the 
elementary elongational instability of a circular flux do- 
main, as illustrated in Fig. |l4|. Again the area fraction 
evolves toward p n — h, but now the deviation is signif- 
icant due to larger contributions from the Biot-Savart 
integral. This shape evolution shows one the means by 
which finite stripes may form in the intermediate state. 




FIG. 15. Elongational instability of a circular flux domain, 
with aspect ratio as in Fig. Il3l but h = 0.38. 
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The curious feature of bulbous tips to the stripe is 
a rather common observation in dipolar systems. It is 
suggestive that the instability is in some sense related 
to evolution toward the fissioning of the original circle 
into two smaller ones. Energetic arguments based on 
this picture show that it rather accurately predicts the 
onset of this instability [p7j . 

The third phenomenon of interest is the branching in- 
stability that occurs for still higher values of h, as shown 
in Fig. |l5|. The initial condition for this simulation was 
a circle perturbed with a small amplitude mode of az- 
imuthal number 3. Rapid growth of that mode is fol- 
lowed by relaxation to "arms" of rather uniform width. 
The angles of the three "arms" forming the vertex are 
close to 120°, as is typical in systems governed by sur- 
face tension. 

We conclude from these studies that a physical mech- 
anism to produce the branched and fingered stationary 
shapes of flux domains in the intermediate state is the 
mechanical instability illustrated in Figures [L3 and Era. 




FIG. 16. Numerical study of branching instability, with 
h = 0.45. 



Finally, Fig. 17 shows the effects of surface tension 



anisotropy on the branching instability of the same initial 
condition as in Fig. 
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While the time evolution first 
produces a four-fold vertex, it subsequently fissions into 
two three-fold vertices that move away from each other. 
The branches of the pattern have oriented themselves 
with respect to the low-tension directions determined by 
the anisotropy (indicated by arrows). The instability of 
vertices of higher order than three is a common feature 
of dipolar systems. 




FIG. 17. Numerical study of branching instability, with 
h — 0.45 and surface tension anisotropy (e = 0.1, q = 4). 
Arrows indicate easy axes of low surface tension. 



V. PERTURBATIONS AROUND THE LAMINAR 
STATE 

In the absence of any in-plane component to the ap- 
plied magnetic field, flux domains in the intermediate 
state often have the shape of buckled lamina. The wave- 
length is typically larger than the stripe width, as is the 
amplitude of the modulation. The conformal mapping 
algorithm for the laminar state is not generalizable to 
treat such truly three-dimensional structures, and there 
does not appear to have been any stability calculation of 
the laminar state. In the following sections we compute 
both the stability and elastic properties of flux stripes, 
making connections with pattern forming properties in 
other systems. 

A. Energy and stability of a single flux stripe 

Consider first the properties of a single flux stripe as 
described by the current-loop model. If the stripe has 
width w in the x-direction, a length I in the y-direction, 
and the plate spacing is d, then by considering the self- 
and mutual induction of the currents flowing along the 
edges the reduced energy E = E/2(tsnA per unit area is 



E 



1 

a 



Nb 
a 



(5.1) 



where a — w/al and Nb = 2M 2 d/asN is the dimension- 
less magnetic Bond number. The integrals are standard 
and yield the result 
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a ia 



ln(l 



2 ) 



+4a tan 1 a - ln(l + a 2 )] 



(5.2) 



Figure [18] shows the stripe energy as a function of its 
width for various Bond numbers. We see that the min- 
imum of this energy becomes sharper as Nb increases. 
Minimizing E with respect to a (hence with respect to 
the width) at fixed area A yields a relation between a 
and Nb, 



l-^[a 2 ln(l + a- 2 ) + ln(l + a 2 )] = 



(5.3) 



Now we connect this result to the stability analysis of 
the stripe. As shown in Fig. |l^(a) and (b) there are two 
classes of small distortions we must consider. The first 
( "peristaltic" ) involves antisymmetric perturbations and 
changes the local stripe width. This will be of higher 
energy than the symmetric (or "serpentine") distortions 
of Fig. |l9|(b) which preserve the width. It is most con- 
venient to calculate the linearized force acting on t he in - 
terface, for which we use the result quoted in Eq. ( 4.11 ), 
written out explicitly, 



-yJC + 2M : 




|r(s) - r(s')| 2 



(5.4) 
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FIG. 18. Stripe energy density as a function of stripe 
width, for various magnetic Bond numbers. 

The details of this perturbation analysis are given in 
Appendix A. After considerable algebra we obtain the 
force associated with a monochromatic perturbation of 



reduced wavevector q = dk; for the serpentine perturba- 
tion, 



aq 



2yl + a 2 



Ko(q) 



F s {q) = q 2 -2A^{ 7B +lll 

+K (aq) - K (Vl + a 2 <z)}, (5.5) 
and for the peristaltic perturbation 

F p (q) = q 2 - 2N B {- 1E + In f^=) ~ M?/2) 

-K (q) - K (aq) + K Q (V 1 + a 2 ?) }. (5.6) 

Let us now look at the limit of small q for serpentine 
perturbations, 

N 



F s (l) = <! 1 - '-f [a 2 



Nb 
64 



ln(l + a- 2 )+ln(l + a 2 )]| (? 2 

|(l + a 2 ) 2 ln(l + a 2 )-a 4 ln(a 2 ) 
-6a 2 (l-^ E )}q 4 

(5.7) 
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FIG. 19. Peristaltic (a) and serpentine (b) perturbations 
of a flux stripe, (c) Illustration of the change in stripe width 
upon a uniform rotation. 

Serendipitously, the condition of stripe equilibrium is 
precisely that which sets the coefficient of q 2 to zero. This 
can be interpreted as a consequence of rotational invari- 
ance. Note first that by assigning the same function 
to the bottom and top edges of the stripe we have main- 
tained the stripe width at w to linear order in £, but not 
at quadratic order. As shown in Fig. U% for a uniform 
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tilt of the layer boundaries, Q v — constant, the width of 

the rotated stripe is w' = w/^Jl + Q — w—(w/2)( 2 + - ■ ■. 

Such a unifor m til t will cost energy through the "bulk" 
term E(w) in 0) as E{w') - E(w) ~ -{w/2)E'(w)g 
where E' = dE/dw. Now, the coefficient of q 2 in (5.7) 
has the interpretation of an effective line tension, associ- 
ated with an energy 



E = KTeff / dyC 2 



(5.8) 



The rotational invariance argument thus shows that the 
apparent surface tension vanishes at the equilibrium 
stripe width. The surviving terms at 0(q 4 ) look like 
bending energy of a rod, 



E 



dyC 



i 



(5.9) 



but this interpretation is spoiled by the term q ln(g), 
whose presence reflects the fundamental nonlocality of 
the magnetic interactions. 
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FIG. 20. Critical Bond number for serpentine instability as 
a function of stripe size. Typical values of w/d are indicated 
for three experimental systems. 

In the case of peristaltic perturbations, the force has 
a finite value as q — ► 0, reflecting the compressibility of 
the stripe. The small-g expansion is 



Fn 



-2N B In 



TT^) + { 



i 



■ [a 2 In a 



-(1 + a 2 ) In Vl + a 2 ] )q 2 + O (q\q A Inq) . (5.10) 

The equilibrium value of the stripe width as a function 
of Bond number, deduced in (5.3), defines the boundary 



of stability of stripes to serpentine perturbations. Figure 
p(i| displays this critical Bond number N B as a function 
of a = w/d. At fixed a, instability occurs with increas- 
ing Ng, and likewise at fixed Ng instability sets in with 
increasing a. In the figure we have illustrated the aspect 
ratios a corresponding not only to type-I superconduc- 
tors but also for typical experiments on magnetic fluids in 
Hele-Shaw flow (with slab thicknesses and stripe widths 
in the millimeter to centimeter range), and for Langmuir 
monolayers (with domains up to tens of microns across 
and a molecular thickness to the layer). 

At its equilibrium width, the energy of small distor- 
tions is positive, vanishing as q — ► 0. Thus is would seem 
not possible to find a field at which the stripe would be 
unstable to a finite-wavelength mode. This conclusion 
assumes that at any applied field the stripe width has 
its equilibrium value (5.3). Fortuitously, the elegant ex- 
perimental observations on buckling instabilities in Lang- 
muir monolayers |^8| have shown us what happens when 
this equilibrium is not reached. Those observations con- 
cerned the dynamics of buckling when the temperature 
was slowly increased. Since these systems are near a crit- 
ical point of phase separation, relatively small changes in 
temperature produce large changes in the density differ- 
ence between the phases (thus altering the discontinuity 
in dipole density A/x) and in the line tension. These, 
of course, directly affect the stability of stripes, quanti- 
fied by the associated electric Bond number. It was ob- 
served that slow temperature ramps produced no buck- 
ling, while rapid heating showed buckling. This suggests 
p9[ that the dependence of stability on ramp rate is as- 
sociated with a competition with mass transport as the 
stripe width adjusts to keep up with the temperature. 
Under rapid ramps, the width is out of equilibrium, yield- 
ing a nonzero (and potentially destabilizing) coefficient of 
q 2 . Turning to the laminar state, this suggests that in the 
early stages of flux penetration such a mismatch between 
the actual and equilibrium widths allows the buckling in- 
stability to occur. 



B. Elastic properties of the laminar state 

By using the CL model we can also examine the elastic 
properties of the laminar state. This is done by consid- 
ering displacements itj (y) of the SN interfaces away from 
the equilibrium laminar phase, as shown in Fig. In 
the long wavelength (continuum) limit, Ui(y) becomes 
the displacement field u(x,y) 7 and the effective elastic 
free energy becomes 



d 2 r 



B 
2 ' 



1 2 

2 U v 



Ki 2 

2 u vv 



(5.11) 



with B the bulk (compressional) modulus and K\ the 
bending modulus. This result applies to serpentine per- 
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turbations of the lamina; the peristaltic perturbations are 
gapped (like optical phonons), and do not contribute to 
the long wavelength properties. The general form of the 
free energy could have been anticipated from the single 
stripe calculations of the previous section; in particular, 
we see that the distortions in the y-direction appear as 
Uy y (or fc 4 |u(k)| 2 in Fourier space), again signifying that 
the effective surface tension is zero. The nonlinear terms 
are required to preserve the rotational invariance of the 
free energy. The free energy, Eq. (5.11), is identical to 



the elastic free energy of a two dimensional smectic liquid 
crystal (mJ. This analogy is quite useful, as the proper- 
ties of two dimensional smectics have been well studied; 
problems such as mechanical instabilities, thermal fluctu- 
ation effects, and defect structures have been considered. 
We expect many of these same phenomena to occur in 
the laminar state. 



u;(yk 



i-l 




FIG. 21. Schematic diagram showing the serpentine per- 
turbations of the domain walls. The displacement of the i th 
domain wall is Ui(y); in the continuum limit this will become 
the displacement field u{x,y). 
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The bulk modulus is 
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In the thick film limit this becomes 
AasNd 

D = . 



(5.14) 



(5.15) 



The bending and bulk moduli may be combined to form 
the length scale A = ^/K\/B, which is a persistence 
length for the distortion of the laminar structure (not 
to be confused with the superconducting penetration 
depth). For thick films this length becomes 



A 2 /a 2 



3 J2m=i sin 2 (mnh)/m 5 
32^ 2 £^ =1 sin 2 (m^/m 3 ' 



(5.16) 



C. Dislocations in the laminar state 

In many of the images of the laminar state [Q one often 
observes edge dislocations, where half of a normal lamina 
has been inserted into the laminar structure (see Fig. |2^ 
for a schematic diagram). 



The elastic moduli for striped phases in two dimen- 
sional ferromagnets with dipolar interactions have been 



calculated in Refs 
rofluids in Ref. 15; 



1|J52|, and for striped phases in fer- 
541]. The calculation for the laminar 
phase is identical, and we include here only the final re- 
sults. The bending modulus is 
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(5.12) 



where the magnetization is M - 
librium spa cing a = yj Ag?//cl 
by Eq. (4.9). In the thick film limit this becomes 



—H c /4ir, and the equi- 
with fch(h,d,a) given 



FIG. 22. A dislocation in the laminar state with a Burger's 
vector of 1. 

By using our elastic theory, we can determine the dis- 
placement field produced by such a dislocation; this prob- 
lem has been studied in the context of 2D smectics |5j| 
and our calculation closely follows Ref. |55|]. We begin 
with the linearized Euler-Lagrange equations for the de- 
fect displacement field u D (x,y), 



ma6'(x)6(y), 



(5.17) 
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where a source term has been added to the right hand 
side to account for the presence of the dislocation, in such 
a way that the line integral of Vu around the dislocation 
is ma (the Burger's vector), with m the number of half 
sheets inserted and a the lamina spacing. Equation (5.17) 
may be solved using Fourier transforms p5[ , with the 
result that 



u fay) = — sgn(x) 




(5.18) 



For an collection of dislocations centered at {r^} with 
strengths {rrii}, we may introduce the dislocation density 



m(r) = y^mj<S(r- Tj), 



(5.19) 



so that the displacement field is obtained by linear su- 
perposition, 



u°(y) = / dV m(r')G(r - r'), 



(5.20) 



with 



G(r — r') = -sgn(x — x') 



erf 



y-y 



/ 4X\x-x'\ / 



(5.21) 

The effective free energy for the d efects, Try, is then 
obtained by substituting Eq. ( 5.20| ) into the the elastic 
free energy, and using the "harmonic conjugate" trick of 
Toner and Nelson fl55||. The final result is 



d 2 r 2 ra(ri)m(r 2 )[/(ri - r 2 ) 



|ri-r 2 |>a 



-E D J d 2 rm 2 (r), 



where the interaction potential is 
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(5.23) 



and the defect core energy is 
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X 2 (GyyY + (G^ 
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(5.24) 



The core energy can be calculated as a function of the 
reduced field h by using the results of the previous section 
[Eqs. (5.14) and (5.16)]. A rough estimate shows that 
this energy is generally of order 10 _3 (_ff 2 /87r)a 3 , and can 
therefore be quite small; as a result it should be easy to 
nucleate dislocations in the laminar phase. 



VI. DISCUSSION 

The free-boundary approach to the energetics and dy- 
namics of the intermediate state has led us to a clearer 
understanding of the shape instabilities of individual flux 
domains and ordered arrays. At the same time, the corre- 
spondence with smectic liquid crystals suggests that phe- 
nomena observed there should have an analog in type-I 
superconductors. Below we suggest several experiments 
to visualize these phenomena. These have as their start- 
ing point the ordered laminar state, produced with an in- 
plane component to the field. Preliminary experiments 
by Reisin and Lipson |56| have shown some of these phe- 
nomena. 

(i) The buckling instability: Rapid changes in the mag- 
nitude of the applied normal field may allow buckling 
instabilities to occur in much the same way as observed 
in Langmuir monolayers. Of interest would be the de- 
pendence of buckling wavelength on the magnitude of 
the field jump. 

(ii) The chevron instability: If the in-plane magnetic field 
is applied at an angle with respect to the lamina we ex- 
pect an instability toward a zig-zag or chevron pattern as 
the lamina attempt to reorient. This is the analog of the 
Helfrich-Hurault effect in smectics f50fl , wherein a field 
component normal to smectic layers, producing a torque 
on them, induces an undulatory instability. 

(iii) The Eckhaus instability: If the magnetic field normal 
to the slab is slowly increased or decreased in magnitude 
the stripe width and spacing must adjust to stay in equi- 
librium. As in convective systems J57|, this may occur 
through an Eckhaus-like nucleation phenomenon to cre- 
ate or destroy laminae. Dislocations can be produced 
which will move toward the sample edges or annihilate 
at the center in accord with the direction in which the 
wavelength must adjust. Their climb and glide dynamics 
will provide an important testing ground for the theory 
of superconductor-normal interface dynamics. 

(iv) Critical-point effects: In the simplified contour dy- 
namics in which the local field is taken to be the critical 
field we saw that the effective Bond number depended 
only on the ratio d/A(T). Near the zero-field critical 
temperature for the superconductor-metal transition the 
interfacial width A(T) diverges with reduced tempera- 
ture (T c — T)/T c . This should produce characteristic 
changes in the equilibrium stripe width as well as pos- 
sibly inducing shape transformations. 

We close by emphasizing what has not been accom- 
plished in this study. First, we have considered laterally 
infinite samples, so that the whole issue of flux penetra- 
tion at the edges is ignored. This is known to be very sig- 
nificant in both type-I and type-II superconductors fl58| . 
A treatment of these effects requires not only the elec- 
tromagnetics of the fields in the neighborhood of the slab 
edges but also consideration of processes such as domain 
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fission. Second, we have presented a highly oversimpli- 
fied dynamical picture in which diffusional instabilities 
are absent and only mechanical ones appear. The inter- 
play between the Mullins-Sekerka and these mechanical 
instabilities appears not to have been considered theo- 
retically and may shed some light on various problems in 
flux domain pattern formation. Third, a large-scale nu- 
merical study of the many-domain problem has not been 
attempted, precluding a clear understanding of the true 
"energy landscape" of this strongly interacting system. 
Fourth, the effects of in-plane components to the applied 
magnetic fields has not been incorporated into the free- 
boundary approach in any quantitative way. This will be 
important for a quantitative understanding of the insta- 
bilities described above. Fifth, the extension of matched 
asymptotic methods used in purely two-dimensional sys- 
tems to the slab geometry in which the intermediate state 
appears has not been developed. A detailed study of this 
point would greatly clarify the free-boundary approach to 
flux domain shapes. Finally, coarse-graining approaches 
to domain dynamics analogous to Otto and Kohn's re- 
cent study of magnetic fluid pattern formation may prove 
quite fruitful 
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APPENDIX A: FLUX-STRIPE STABILITY 
CALCULATION 

Here we collect some details of the stability analysis for 
single flux stripes. In the case of peristaltic distortions, 



r±(y) = ye y + [±— ±C(y) j e x , 



(Al) 



The linearized normal force can be reduced to the form 
U (y) = jCyy 

poo 

+2M 2 / dy' [C(y') + C(y)] [Sr+rf - S w ,\ 



+2M 2 / dy' {((y') - av)] [Sd* - S ] , (A2) 



where S a 2 = l/\(y' — y) 2 +a 2 \. For serpentine distortions, 
the displacements are 



r± (y) = ye y + (±^ + C(y)J e x , 

and the normal force has the simpler form 
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If the distortion is the plane wave £(?/) = Acos(ky), then 
the force is proportional to £ with a coefficient F(k) 



F(k) 



- 7 fc 2 + 



4M 2 / dy [1 - cos (ky)} 
Jo 

x [Sd 2 +w 2 — S w 2 — S d 2 + S ] (A5) 



where now S a 2 = l/\y 2 + a 2 \. Several of the integrals 
with integrands proportional to cos(ky) reduce trivially 
to Bessel functions, but care must be taken to account 
for canceling divergences in the remaining terms. This is 
conveniently done by considering the limiting process 
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Rescaling these equations, we obtain 
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A useful intermediate result at this stage is |6l[] 

cos(y) 



1 - cos(y) 
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(A7) 



(A8) 



where 7^ = 0.577215 ... is Euler's constant. Substitution 
into ( A7) and evaluation of the remaining integrals then 
yields the final results ( [3.5] ) and (5.6) for the energy of 
serpentine and peristaltic perturbations. 
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TABLE I. Analogies between free streamline flow in fluids and lamina formation in the intermediate state of type-I super- 
conductors. 



Free streamline flow around a plate 



Complex potential w = 4> + itp 

Complex fluid velocity u — iv = 

Streamlines 

Free streamline 

Free streamline velocity U 

Region of fluid flow 

Cavity behind plate 



-dw/d( 



Rial 



flc 



Laminae in superconductors 



Complex potential w = <f> + iA y 

Complex magnetic field B — B x — iB z — - 

Field lines (lines of force) 

Superconducting-normal interface 

Superconducting critical field H c 

Normal phase with nonzero magnetic field 

Superconducting phase 

Lamina in a finite thickness plate 



dw/d( 
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TABLE II. Analogies between interfacial energetics of type-I superconductors and other systems. The energy of a set {T>i} 
of domains is written as _E[{i"i}] = Ai + r y^2i i Li — <j>ds <j>ds'ii ■ tj$y {Rij/Q ■ 



System 


n 


7 




$ 




£. 


type-I superconductors 


(HZd/8ir){ Pn + h*/p n ) 


H 2 dA/8n 


H 2 d/8iv 2 
2dM 2 


sinh 




d 


magnetic fluids 6 


Lagrange multiplier 


d(J F\y 


sinh -1 


'l/z)+z-Vl + z 2 


d 


Langmuir monolayers 2 
FitzHugh-Nagumo modeL 


Lagrange multiplier 
AF 


JLE-LC 

D 


(Am) 2 
P 


1/22* 

K (z) 




l" 



Explanation of symbols: o~FWi ferrofiuid water surface tension; M, ferrofluid magnetization; Jle-lc, line tension 
between liquid expanded (LE) and liquid condensed (LC) phases in a Langmuir monolayer; A/i, discontinuity in 
electric dipole moment density between LE and LC phases; d mo \, a molecular cutoff - monolayer thickness. 
a Present work. 
6 Refs. 

c Refs. J||7|. 

d Eq. Q and Refs. H|o|. 

* This limiting form is supple ment ed with a cutoff procedure. See Ref. [B7| . 
** The system of units in Eq. (|]^) sets the inhibitor screening length to unity. 
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